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PROPAGATION AND STABILITY OF SUPERLUMINAL WAVES IN PULSAR WINDS 
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ABSTRACT 

Nonlinear electromagnetic waves with superluminal phase velocity can propagate in the winds around 
isolated pulsars, and around some pulsars in binary systems. Using a short-wavelength approximation, 
we find and analyze an intcgrable system of equations that govern their evolution in spherical geometry. 
A confined mode is identified that stagnates to finite pressure at large radius and can form a precursor 
to the termination shock. Using a simplified criterion, we find this mode is stable for most isolated 
pulsars, but may be unstable if the external pressure is high, such as in the pulsar wind nebulae in 
starburst galaxies and in W44. Pulsar winds in eccentric binary systems, such as PSR 1259-63, may 
go through phases with stable and unstable electromagnetic precursors, as well as phases in which the 
density is too high for these modes to propagate. 

Keywords: plasmas - instabilities - waves - pulsars: general - stars: winds, outflows - ISM: supernova 
remnants 



1. INTRODUCTION 

Electromagnetic fields, modulated at the rotation fre- 
quency of the neutron star, form the energetically dom- 
inant component of pulsar winds. These flows are re- 
sponsible for transporting the rotational energy lost by 
the star and depositing it in the surrounding pulsar wind 
nebula (PWN). As well as energy, they also convey the 
magnetic flux and the charged particles — most likely 
electrons and positrons — that are required to produc e 
synchrotron radiation in the PWN (jRees fc Gunnlll97l ). 

An important property that sets pulsar winds apart 
from other stellar winds is their relatively low density. 
As a result, the fluctuations imposed by the rotation of 
the neutron star are able to propagate not only as MHD 
waves frozen in to the outflowing plasma, but, beyond 
a critical or cut-off radius r c , also as large-amplitude 
electromagnetic modes of superluminal ph ase velocity 
(|Asseo et alj|1975t iMelatos fc MelrosdfT996l) . The loca- 
tion of the cut-off radius depends on the wave ampli- 
tude and on the relative strength of the phase-averaged 
fields and the fluctuating components, which, in turn, 
depend on the obliquity of the pulsar (i.e., the angle be- 
tween its magnetic and rotation axes) and on latitude 
in the wind. For isolated pulsars it generally lies well 
inside t he position where a termination shock can be ex- 
pected (jArka fc KirkI 120 12t ). so that superluminal waves 
may play an important role in the outer parts of the 
wind and also in the termination shock itself. In particu- 
lar, since they may carry the entire wind luminosity and 
are notorio usly unstable (lMaxlll973t [Drake et all 1 19741 : 
lAsseo et all 119801: iLee fc Lerchel 119801) . thev could pro- 
vide the key to resolving the well-known "er-problem" . 

This term is used to describe the lack of a convinc- 
ing mechanism for converting Poynting flux into particle- 
carried energy flux. Poynting flux dominates the wind at 
launch, but is thought to be a small fraction of the en- 
ergy budget outsid e the termination shock. Recent work 
(jPorth et al.l2012t ) suggests that M HD instabil i ties in the 
shocked wind may, as proposed by IBegelmanl (|1998l ) , be 
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able to reduce a (the ratio of Poynting flux to kinetic- 
energy flux) from a starting value of a few just outside 
the shock, to the value sugges ted by observation, which i s 
- 10~ 3 for the Crab Nebula (jKennel fc CoronitillT984rl . 

However, the mechanism underlying the transition 
from a 3> 1 near the pulsar to a ~ 1 at the ter- 
mination shock remains controversial. Most work on 
this problem has concentrated on th e damping of MHD 
like wave-modes either in the wind dLyubarskv fc KirkI 
[2001 IKirk fc Ski n1 l200l ILvubarskvl 1201(1 " or at 

the termination shock itself, where current sheets car- 
ried into the shock are compressed, giving rise t o dis- 
sipation by driven reconnection. |Lvubarskv| | 2003t 
iLvubarskv fc Livertsll2008HPetri fc Lvubarskvll2007t) . In 
this scenario, particle acceleration is primarily asso- 
ciated with the dynamics of the reconnection region 
( Sironi fc Smtkovsk^l2llTl . 

However, it has recently emerged that the termina- 
tion shock structure may be very different in the region 
where superluminal modes can propagate. In this case, 
the mechanism of dissipation is closely connected with 
the parametric instabiliti es of these modes, ra ther than 
with driven reconnection (jAmano fc Kirldl20r! . The as- 
sociated particle acceleration mechanisms have not yet 
been investigated in detail, but it is suggested that the 
first order Fermi process may be much more prominent 
when superluminal waves are present, than in the case of 
dissipation by reconnection. 

This question is central to the study of PWN, and es- 
pecially their high energy (TeV) emission. Current mod- 
els adopt ad hoc assumptions concerning the injection 
of accelerated particles at the termination shock which 
take no account of the role of the surrounding medium 
in determining the shock location and structure. In or- 
der to establish a predictive theory of pulsar winds, it 
is, therefore, fundamentally important to identify which 
PWN are likely to sustain shocks mediated by superlu- 
minal modes, and which are not, and to classify impor- 
tant properties such as the stability or instability of these 
modes at the position of the shock. These are the ques- 
tions we address. Our approach builds on the work of 
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lArka fe Kirkl (|2012f ). who found the shape of the cut-off 
surface r — r c for both linearly and circularly polarized 
waves of arbitrary amplitude, including, for the former, 
a non-zero value of the phase-averaged magnetic field. 
Here, we develop a perturbation method to determine 
how a superluminal mode launched outside the cut-off 
surface will evolve as it propagates radially outwards, 
searching for those regions of the wind in which such a 
packet remains relatively stable, and those in which it 
can be expected to thermalize rapidly. 

In section [2] we recall the description of nonlinear su- 
perluminal plane waves in a two-fluid electron-positron 
plasma, and summarize the literature concerning their 
stability properties. Radially propagating waves in 
spherical geometry are treated in in section |3l First, the 
short-wavelength perturbation theory is developed and 
is shown to result in an integrable system of equations 
for the radial evolution of wave packets, consisting of the 
conservation laws of particles and energy, supplemented 
by an analogue of the entropy equation. Restricting, for 
simplicity, the treatment to circularly polarized modes, 
it is shown that two kinds of wavepacket are possible — 
freely expanding modes, which ultimately turn into vac- 
uum waves, and confined modes, which stagnate at finite 
pressure, and can be identified as extended precursors of 
the termination shock. Section 13.41 discusses the con- 
fined modes and shows that two distinct regions exist in 
which they are relatively stable. One of these is relevant 
for pulsar winds in high pressure environments, such as 
the high-pressure wind of a companion star, the other 
is relevant for isolated pulsars in low-pressure surround- 
ings, such as supernova ejecta. Analytical constraints are 
given on the values of pulsar period, surrounding pressure 
and mass-loading parameter fj, that determine whether or 
not a pulsar wind has the potential to generate a stable 
superluminal precursor at its termination shock. In sec- 
tion U we discuss the significance and limitations of our 
results and their implications for PWN around isolated 
pulsars and in gamma-ray binaries. Finally, a brief sum- 
mary of our main results is presented in section [5] 

2. LARGE-AMPLITUDE PLANE WAVES 

2.1. Basic properties 

The simplest system able to describe superluminal 
modes is that of t wo cold charged fluids. I t has been ana- 
lyzed in detail bv lClemmowl (|1974al fl977T l . In the pulsar 
case, we are concerned with cold electron and positron 
fluids and, in particular, with electromagnetic modes in 
which the displacement current is non-zero in all frames 
of reference. 

In the outer parts of a pulsar wind, the radial com- 
ponent of the magnetic field is very small, and will be 
neglected in the following. In this case, radially propa- 
gating electromagnetic waves have only transverse fields, 
but the fluid velocities can have a component in the di- 
rection of propagation. The waves are characterized by 
a superluminal phase speed /3 p hasc > 1 and a sublumi- 
nal group speed /3 W = l//3 p hase and will henceforth be 
called simply "superluminal" modes. The electron and 
positron fluids have the same proper number density n, 
and move with the same velocity P11/7 (in units of c) in 
the direction of wave motion. But the transverse velocity 
of the positrons p±/j is antiparallel to that of the elec- 



trons. (Here, 7 is the Lorentz factor of the fluids, and p\\ 
and p± are components of the dimensionless four- velocity 
of the positron fluid). This generates a transverse con- 
duction current. In the homogeneous or "H- frame", in 
which the group speed vanishes, the fluid and field vari- 
ables are space- independent, and the conduction current 
exactly balances the displacement current. In general, 
p\\/"f differs from the wave group speed /3 W , implying a 
non-vanishing flux of particles in the H-frame. 

The g overning equations have been presented else- 
where fin lArka fe K irk 20f 2, for example), and are repro- 
duced in Appendix [XJ The nonlinear solutions take their 
simplest form in the H-frame, where the phase variable 
(|A1|) is purely temporal. Circularly polarized modes are 
monochromatic (i.e., depend sinusoidally on phase) and 
have phase- independent pn. The phase-averaged fields 
vanish. Linearly polarized modes, on the other hand, 
have a "saw-tooth" dependence of on phase, and may 
also carry a non-zero phase-averaged transverse magnetic 
field component perpendicular to the oscillating electric 
field. 

Although the analytical treatment of linearly polarized 
modes is significantly more cumbersome, their dispersion 
curves, at least for vanishing phase-averaged magnetic 
field, a re very similar to those of circularly polarized 
modes (jArka fc Kirkl I2012T ). Consequently, we restrict 
much of the following discussion to circularly polarized 
modes. 

2.2. Stability 

In unmagnetized plasmas, strong electromagnetic 
waves are subject to parametric instabilities. These 
are induced by longitudinal density fluctuations which 
can couple to the transverse electromagnetic side-band 
modes and cause ba ckscattering, filamentation, or ab- 
sorption of the driver (|Maxlll973l : lDrake et al.lll974l ). The 
growth rates can be as high as the wave frequency. The 
stability of self-consistent waves was fi rst investigated 
in a nonrelativis tic electron-ion plasma ([Max fc Perkins! 
ElZllMllIIiZi)- This work was extended to the rela- 
tivistic cas e for electron-positro n plasmas by IRomeirasI 
(|1978f l and iLee fc Lerchd ()1978f ). In all cases only per- 
turbations that propagate in the direction of driver's mo- 
tion were discussed, and the dispersion relations were 
obtained by linearization. In general, it was found that 
both short- and long-wavelength perturbations are un- 
stable. Finite-temperature effects can be expected to 
suppress the instabilities at short wavelengths, although, 
to our knowledge, a complete analysis is lacking. Long 
wavelength perturbations, on the other hand, are stable if 
the particle flux through the wave in the H -frame is suffi- 
ciently large. A simple criterion is given bv lLee fc Lerchel 
(|1978f l. Denoting quantities measured in the H-frame by 
a prime, they find stability when 

S=p\f-2 1 'p ± +pl > . (1) 

Numeric al studies of IRomeirasI (11978f ) and the PIC simu- 
lations ofrSkia:raa.s cn et al.l (|2005j ) confirm the stabilizing 
effect of relativistic streaming. Although no precise test 
of it has been undertaken, we nevertheless adopt ([T]) in 
the following as the stability criterion for circularly po- 
larized modes, noting that for linearly polarized modes 
the presence of a phase-averaged transverse component 
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of the magnetic fie ld has an additional stabilizing effect 
(jAsseo et alj|1980j) . 

2.3. Electromagnetic Hugoniot curve 

Plane waves in this model can be characterized by the 
phase-averaged values of their particle flux density J, 
energy flux per particle /i, and parallel (i.e., in the prop- 
agation direction) momentum flux per particle v. In a 
local analysis, waves in a radial pulsar wind can also be 
considered plane. A radial wind occupying a total solid 
angle f2 s is characterized by an energy flux L and a flux 
of electrons and positrons of N, in terms of which the 
parameters \i and J are 



(2) 
(3) 



li = L I [Nmc 2 
J = N/ (n s r 2 ) 



The momentum flux is connected with the magnetizati on 
parameter a. In a cold, striped wind, it is (lKirkll20iof ) 



a 

2fl 



8n 3 



(4) 



(for ji ^> 1 and a < /i 2 ^ 3 ) and is independent of radius, 
provided there is no dissipation. It is related to the ram 
pressure P (the (r, r) component of the stress-energy ten- 
sor) by 



P-- 



L 



fj, \Q s r 2 c 



(5) 



Observations of the Crab Nebula suggest (e.g. 
iBucciantini et alJl201l| ) 



L = 5 x 10 JB erg s" 

s 



A>=10 40 " _1 



(6) 
(7) 



so that /i ps 10 4 . Close to the star, a magnetically 
dominated wind is expected to accelerate rap idly until 
mildly super (magneto)sonic (jKirk et al.l I2009T ) implying 
< 2 / 3 so that v ~ ll. For nonlinear waves, the wave 
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group speed depends not only on the frequency, but also 
on other wave properties, such as the amplitude. Thus, 
in order to plot a dispersion curve, giving, for example 
the group speed as a function of frequency, additional 
constraints are needed. Choosing \i and v to be constant 
along such a curve transforms it into the electromagnetic 
equivalent of a Hugoniot curve: it then specifies all the 
(plane) waves of a given frequency that can be launched 
in a p ulsar wind characte rized by the same values of \i 
and v (lArka fc Kirkll2012h . 

3. SPHERICAL WAVES 

3.1. Short-wavelength approximation 

The radial evolution of spherical nonlinear waves 
at distances from the origin large compared to the 
wavelength can be treated us i ng standard perturbat ion 
techniques (jAsseo et all 11984 : iKirk fc Mocho]||2011allbl ). 
These lead straightforwardly to equations that express 
conservation of the phase-averaged particle and radial 



energy fluxes (Appendix [A"]) 

_t__d 

r 2 dr 

[_d_l r2 P- 
r 2 dr 



l {2np V ) 




Airmc 2 



(8) 
(9) 



However, the divergence of the (r, r) component of the 
stress-energy tensor does not yield a conservation law 
directly: 




(npl) 



(10) 



Nevertheless, as we show in the Appendix [B] an integral 
of motion can be found by constructing the electrody- 
namic equivalent of the hydrodynamic entropy equation, 
essentially by subtracting /3 W times the energy equation 
(|9|) from the momentum equation (|10[) . This integral can 
be found more directly by considering the adiabatic in- 
variants of the zeroth-order plane waveQ In both the lab. 
(pulsar) frame and the H-frame, the charge density of the 
wave vanishes, and both the current and the electric field 
are entirely transverse. This means that the electrostatic 
potential A vanishes when the Coulomb gauge is cho- 
sen. It is then a trivial matter to construct an adiabatic 
invariant $ of the motion of an individual fluid particle 
by integrating over one period of oscillation the zcroth 
component of the canonical momentum P^: 



dtP° 



= I dt (mc7 + eA°/c) 



- 2innc (7) /lj 



(11) 



where u) is the angular wave frequency in the lab. frame. 
In our case, ui equals the (constant) angular velocity of 
the pulsar, so that the invariance of $ in a slowly ex- 
panding (or contracting) flow implies 



d 

^-(7>=0 
or 



(12) 



In addition, to these equations, Ampere's law provides 
a connection between the electric field and the fluid mo- 
mentum: 



— = -mcuj p -/ wP± /e , 

-V2 • 



(13) 



where 7 W = (l — ft 2 ,) " is the Lorentz factor associated 

1/2 

with the wave group speed, and u> p = (8irne 2 j 'mj is 
the proper plasma frequency. 

3.2. Circular polarization 

For simplicity, we now restrict the treatment to circu- 
lar polarization, the corresponding expressions for linear 
polarization are given in Appendix [Cl 

We are indebted to an anonymous referee for suggesting this 
approach. 
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To close the system, a non-linear dispersion relation is 
required. For circular polarization, in which the variables 
n, 7, \E\ 2 , |pj_| 2 , andp|| are all independent of phase, this 



2 2 

: 7wW p 



and equation (|T3j) reduces to 
e 2 \E\ 2 



\PA 



(14) 



(15) 



Equations ([5]) and can be integrated and combined 
with dHJ) and {jT5j) to give 



J = r 2 2np\\ 



M = 7 + 



7wftwPj 



(16) 
(17) 



where J and [i are the constants of integration. Equa- 
tion (|17|) simply expresses the conservation of energy, in 
which the first term is the particle energy flux and the 
second the Poynting flux (per particle, in each case). In 
a pulsar wind, both of these quantities, and also the par- 
ticle flux J, are expected to be positive, corresponding 
to outward-going fluxes. Thus, the relevant parameter 
space is restricted to p|| > 0, /3 W > 0, and 7 < /1, and the 
bounding line fi = 7 corresponds to waves of vanishing 
amplitude, p± = 0. 

As a wave packet propagates, its frequency remains 
locked to that of the pulsar. This determines the radial 
dependence of the wave via the continuity equation f| 16[) 
and the dispersion relation (|14[) : p|| oc 7^,/r 2 . In terms of 
the normalized radius R defined in lArka fc Kirkl ((20121 
Eq. 31): 



R = r{tt s Loj 2 /AncN 2 e 2 ) 1/2 , 



this relation is 



M7w 
R 2 



(18) 



(19) 



(Note that the cut-off radius r c lies very close to, but just 
outside, the point R = 1.) 

Using 7 2 = 1 + p 2 + p\ to eliminate p± and equa- 
tion (|19[) to eliminate p|| , we rewrite the equation of con- 
servation of energy flux per particle (fT7|) as 



i? 4 



+ (M-7)^-^(7 2 -l)=0 • (20) 



This equation determines the radial dependence of the 
group speed of the wave /3 W , because, according to the 
entropy equation (fl"2"|) . 7 = constant. The radial depen- 
dence of p±, \E\ and n follow from (HHJ), dTTJ) , (fl"5|) 
and p6|) . The ram pressure, normalized to the value of 



L/n s 



at the cut-off radius is 



P- 



liR 2 



(21) 



where v w (R), which is the ratio of the momentum flux 
density to the particle flux density, is defined in analogy 



with ([T7J) as 



i/ w (R)=p\\ 



ilv\ (1 



2p,| 



(22) 



For circular polarization, Eq. (f2"0")) also allows one to 
construct the electr omagnetic Hugoniot curve found by 
lArka fc Kirkl (|2012f ). Whereas the dispersion curves fol- 
low from the additional constraint 7 = constant, the 
Hugoniot curve is the locus of points under the con- 
straint v w (R) = v = constant. In the limit of a flow 
carrying zero Poynting flux, which can convert only into 
a wave of zero amplitude, these two curves coincide, 
7 = /j, = v w = constant, and lie on the line 



1/2 



(M 2 - 1) 



-1/4 



(23) 



3.3. Confined and freely expanding modes 



The radial evolution of p± and P, for [i = 10 4 arc 
shown in Fig. Q] for various values of 7, together with the 
Hugoniot curve for a particular choice of a (and, hence, 
v). 

At large R, two modes propagate: a freely expanding 
mode in which the second term in (|20|) is negligible: 



pii^(7 2 -i) 1/2 , 
7w ^( 7 2 -i) 1/4 j r/Vm , 



p- 



^-7+y / 7 2 -l 
fiR 2 



(24) 

(25) 
(26) 

(27) 
(28) 



and a confined mode in which the first term in (|20[) is 
negligible: 



P^(7 2 -l) 1/2 

o . Mm -7) 

Pv 



( 7 2 -l)i? 2 ' 
7w^l , 

7 2 -l 



P- 



2/i 2 



(29) 
(30) 

(31) 

(32) 

(33) 



In the free-expansion mode, the direction of the fluid ve- 
locity aligns itself with the radial direction as the mode 
moves outwards. The radial fluid velocity tends to a con- 
stant value, the density and ram pressure drop as 1/i? 2 , 
and the wave amplitude decreases as 1/i?. In the con- 
fined mode, on the other hand, the fluid velocity turns 
towards the transverse plane. The radial particle flow 
stagnates, and the ram pressure, density and wave am- 
plitude all tend to constant values, such that the proper 
plasma frequency equals the pulsar rotation frequency. 

In a wind-nebula system, in which the outflow is con- 
fined by the external medium in a slowly expanding bub- 
ble, the wind must terminate by decelerating to nonrel- 
ativistic speed at a point where its ram pressure roughly 
equals the external pressure. This is achieved at a shock 
discontinuity in an MHD picture. From Fig. [1] and 
Eq. ([33]) . we see that the same effect is produced if the 
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Figure 1. The electromagnetic Hugoniot curve (red) for = 10 4 , 
a = 100 and the radially propagating wave modes (blue) for cir- 
cularly polarized nonlinear waves launched at the blue dots. Top 
panel: p± (lower branch: free-escape mode; higher branch: con- 
fined mode). Bottom panel: the ram pressure P, defined in Eq. i ll' 1 1 
for the confined mode only. 

flow converts into a confined superluminal mode, which 
then decelerates at almost constant ram pressure. If 
this mode remained stable, it could in principle accumu- 
late and fill a large volume around the conversion point, 
which would be see n as a pulsa r wind nebula, as orig- 
inally suggested by iReesI (|1971| ). The location of the 
conversion point is determined by the external pressure: 
the higher the surrounding pressure, the closer it lies to 
the pulsar. As can be seen from Fig. [2] the ram pressure 
changes by less than a factor of 2 between the conversion 
point and R — > oo. Consequently, this point is located at 
roughly the position expected for the termination shock 
in the MHD picture. 

3.4. Zones of stability 

However, as discussed above, once the mode has slowed 
down sufficiently, it is subject to strong parametric in- 
stabilities, which ultimately dissipate the coherent oscil- 
lation, leaving behind the relativistic electron positron 
plasma and residual magnetic field that fill the pulsar 
wind nebula. In analogy with the MHD picture, we refer 
to the point at which instabilities arise in this scenario as 
the "termination shock" . The region between the point 
at which the confined mode is launched and the termi- 
nation shock is then an extended "precursor" . 

The propagation and stability properties of the precur- 
sor are completely determined by equations ([20]) and (TTJ 
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Figure 2. The fractional difference between the ram pressure 
^conv at the conversion point and its value P^ at R — v oo, plotted 
for the confined mode with fi = 10 4 and a = 100. 

and can be extracted analytically and plotted relatively 
straightforwardly. In pulsars, we are concerned with rel- 
ativistic flows of low mass- loading, fi 3> 1, which can 
be assumed to be supersonic a < ^i 2 / 3 . Therefore, we 
present here only simplified expressions obtained to low- 
est order in an expansion in powers of e ~ /^ _1 / 3 ~ er~ 1//2 . 
In this approximation, the stability criterion (fT]) for a sta- 
tionary wave mode (/3 W = 0) is 

S^n 2 (i? 2 - 2V-R 4 - l) /R 2 > . (34) 

Thus, stationary modes, for which the H-framc coincides 
with the laboratory frame, arc stable only close to the 
cut-off radius, in the range 1 < R < (4/3) 1 / 4 . This region 
is illustrated in Fig. [3] In fact, the region of stability 
extends also to modes which are not stationary in the 
lab. frame. As is clear from the definition of S in Eq. (fT]), 
the waves arc formally stable along the limiting line for 
weak waves (|23|) . on which p± = 0. But the region of 
stability reduces to a thin layer in the neighborhood of 
this line as R increases. Assuming 7 W ~ R > y/a, we 
find approximate roots of S at 

i? = 7w (8±4V3) . (35) 

Thus, at large 7 W , the zone of stability lies in the range 
7 W < R < 1.077 w . 

Equation ([33)) reveals a second stable region, which, for 
large 7w , lies at R > 14.9 7w . This region is illustrated in 
Fig. |U For R ~ a, it is bounded from below by the line 
7 W = 2, but this boundary drops to lower wave group 
speeds at larger R. The stable zone vanishes inside a 
critical radius, which is roughly where the line 7w = 2 
intersects the line R = 7 W (8 + 4\/3) , at R = 30. 

The Hugoniot curves for a = 100 are also shown in 
Figs. [3] and 2) Expanding the expression for the con- 
fined mode branch to lowest order in e ~ /i" 1 / 3 ~ er -1 / 2 , 
and combining this with the approximate expression for 
the lower bound of the stability curve (see Appendix [Dj 
shows that the launched waves are stable, provided the 
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Figure 3. The inner zone of stability. Electromagnetic Hugo- 
niot curves arc plotted showing the group four-speed (/3 W 7 W ) as a 
function of radius R for fi = 10 4 , a = 100 (red lines, solid: free 
expansion branch, dashed: confined branch). In blue, the radial 
evolution of three confined modes launched at different radii are 
plotted. In the green region the waves arc stable according to the 
criterion given in Eq. Jl}. In the pink shaded region they are un- 
stable. The unshaded region depicts waves with inwardly directed 
Poynting flux which are not relevant to pulsar winds. 

launching radius exceeds a critical value 

R > -Rcrit 

wlOO . (36) 

The length of a precursor consisting of a stable, nonlin- 
ear electromagnetic wave can be estimated from Figs [3] 
and [H In the inner stable region, the slowing down of 
the wave plays a relatively minor role. Even at con- 
stant group speed, the precursor wave enters the unsta- 
ble zone after propagating for at most 4% of the radius. 
Depending on the pulsar parameters, which fix the re- 
lationship of the dimensionlcss scaled radius R to the 
light-cylinder radius, the stable precursor can neverthe- 
less extend over many wavelengths. In the outer stable 
region, however, it is the wave deceleration that drives 
the precursor wave into the unstable region. The ex- 
tent in radius over which this occurs can be estimated 
as AR/R w dlogi?/dlog(/3 w 7 w ). Using Eq. {20]), we 
find for the confined mode AR/R w 1/8, so that here 
again, the precursor extends over a relatively small frac- 
tion (~ 10%) of the radius, which may, however, corre- 
spond to many wavelengths. 

4. DISCUSSION 

The effects of spherical geometry on strong, radially 
propagating electromagnetic modes are potentially im- 
portant for the structure of pulsar winds and their ter- 
mination shocks. However, pr evious work on this topic 
(|Asseo et all 119841: lKirkll2010l ). is flawed because it as- 
sumes the divergence of the radial momentum flux van- 
ishes. In fact, as we show above and in Appendix [Bj the 
radial momentum equation does not yield a conservation 
equation; as the modes propagate, the transverse and 
parallel degrees of freedom exchange energy with each 
other, keeping the total energy is conserved. We show 
that, for arbitrary polarization, there exists an additional 
conserved quantity — the phase-averaged Lorentz factor 




12 3 4 

log(R) 



Figure 4. The outer zone of stability. The same Hugoniot curves 
are plotted as in Fig. [3] together with three confined modes, 
launched at larger radius. The pink and green shadings again de- 
pict unstable and stable regions, respectively. The inner zone of 
stability lies close to the Hugoniot curve of the free escape mode, 
and is not visible on the scale of this figure. 

of the particles (7) — which, together with the conser- 
vation of energy and particle fluxes, enables the system 
to be integrated. 

Nevertheless, the case of general polarization remains 
unwieldy, so that the discussion we present of the mode 
properties is mainly restricted to circularly polarized 
modes, for which phase-averages can be dropped. The 
dispersion curves, which we h ere call electromagn etic 
Hugoniot curves, were found bv lArka fc Kirkl ([2013 ) for 
linearly polarized modes with vanishing phase-averaged 
magnetic field. They closely resemble those for circular 
polarization, suggesting that this simplification is not an 
important restriction, at least in the region close to the 
equatorial plane. 

The electromagnetic Hugo niot curve admits t wo solu- 
tions for radial propagation (jArka fc Kirld[20T2T ) . By ex- 
amining the asymptotic behavior of waves launched on 
these branches, we identify them as a free-escape and a 
confined mode, with very different properties. The group 
speed of the former accelerates, whereas that of the latter 
decelerates during outward propagation from the launch- 
ing point. We argue that, in outflows such as those from 
pulsars surrounded by nebulae, only the confined mode 
can provide a self-consistent solution matching a rela- 
tivistic, highly magnetized wind to a slowly expanding, 
weakly magnetized nebula. We show that the external 
pressure determines the radius at which a confined mode 
must be launched in order to realize such a solution. The 
higher the pressure, the closer to the star the launching, 
or conversion, takes place. However, we do not discuss 
the physical mechanism by which the wind converts from 
one mode to the other. This question demands a different 
approach, and most likely necessitates either two-fluid 
or part icle-in-cell simulation s. It has been addressed re- 
cently (jAmano fc K irk 2013), but only for a very limited 
range of physical parameters, corresponding to a launch- 
ing radius close to the critical radius r c . 

A confined mode becomes unstable against longitu- 
dinal density perturbations when the particle stream- 
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Table 1 

Pulsar wind nebulae studied by Bucciantini ct al. (2011). 



References 



pulsar 

,38 



(ms) 



L(xl0 38 erg s _1 ) 
a L (xl0 10 ) 
ft (lower limit) 
(j, (upper limit) 
Pconv (dyn cm -2 ) 
Conversion radius 
-Rconv (upper limit) 



Crab 


3C 58 


B1509-58 


Kes75 


W44 


K2/K3 Kookaburra 


1 


2, 3 


4,5 


6 


7, 8 


9 


33 


65 


150 


324.8 


267 


68.2 


5 


0.27 


0.18 


0.083 


0.0043 


0.1 


7.6 


1.8 


1.4 


0.98 


0.22 


1.1 


10 6 


5 x 10 5 


3 x 10 5 


10 5 


10 5 


10 5 


1.9 x 10 4 


8.8 x 10 3 


1.2 x 10 4 


2.4 x 10 4 


5.5 x 10 3 


2.7 x 10 4 


~ 10- 8 


3.2 x lO" 10 


~ 10" 10 


~ ur 11 


> 5.7 x IO" 10 


4 x 10" 10 


575 


781 


800 


2300 


86 


1979 



References . - (1) Kennel k Coroniti 198 4al; (2) [Mu rray ct al. 2001; (3) ISlane et al.l [20021: (4) 

Gacu sler et "all [19991 ; (5) lYatsu et al.ll2009t (6) INg et al.ll2008l ; (7) IFrail et al.N199l (8) ICox et al.M1999k (9) 
Ng et al.1120051 



ing through the wave drops below a critical value. For 
plane, circularly polarized waves, we adopt a simple 
stability condition (JlJ ag ainst long-wavelengt h, parallel 
fluctuations proposed by ILee k, Lerc hc (1978)). In fact, 
obliquely propagating fluctuations may change this sim- 
ple condition if th ey turn out to be m ore unstable than 
the parallel ones (|Lee fe Lerchel [i"978() . but we are not 
aware of a suitable alternative criterion. On the other 
hand, waves which have a nonvanishing component of the 
transverse magnetic field, i.e., those which are launched 
at higher latitudes in puls ar winds, may be stabilized by 
the presence of this field {A ssco "etanfl98(il) . 

Using this criterion, we find two regions in which con- 
fined modes are stable when launched, although they 
subsequently propagate into an unstable region. In 
terms of the normalized radius, these arc an outer zone, 
R > 100, and an inner zone R ss 1. These zones and 
the location of their boundaries are essentially indepen- 
dent of the parameters /i and a of the pulsar wind, pro- 
vided it is relativistic, // > 1, and supermagnetosonic, 



a < /i 



2/3 



However, their relevance depends on where 



the pressure confining a particular pulsar wind places 
the launching point. This, in turn, depends on the ra- 
dius normalization, which varies from pulsar to pulsar. 
In terms of the light-cylinder radius rx, the normalized 
radius is R = (r/r^di/a-^ rj r-^/r c ), where the strength 
parameter at the light cylinder ol is a dimensionless 
measure of the spin-down power per unit solid angle: 

a L = 3.4 x 10 10 [(i/lO^ergs" 1 ) (4rr/fi,)] 1/2 . 

We show in Fig. [2] that the pressure at the conversion 
point lies close to L/fl s r 2 c, which can be written as 



p — . 

1 ronv 



io- 



R? 



10 4 



P, 



nd) 



pulsar 



dyn cm 2 ,(37) 



where P pu isar is the rotation period of the pulsar. This 
must be compared with the nebular pressure, usually es- 
timated from the volume-averaged, equipartition mag- 
netic and particle energy densities r equired to produce 
the n ebular synchrotron emission ([Kennel fe Coronitil 
Il984al) . These estimates, together with the implied nor- 
malized radius of the conversion point R CO nv are listed 
in Table [J for those p ulsar wind nebulae studied by 
iBucciantini et all ([201 1[ ) . We see from this table that the 
isolated pulsars can be expected to launch stable precur- 
sor waves at the termination shock, since in all but one 
case i?conv > 100. The only exception is the pulsar in 



the supernova remnant W44, which, however , is known 
to be interacting with a dense molecular cloud (|Cox et al.l 
[1999ft . 

The PWNe of isolated pulsars in starburst galaxies 
have recently been proposed as sources of the very high 
energy emission from thos e galaxies ([Mannheim et all 
120121 lOhm fc Hinton|[20l3[ ). The pressure in the inter- 
stellar medium in these objects (~ 10 _9 dyncm~ 2 ) is 
some three orders of magnitude greater than in the Milky 
Way. If this scaling is reflected in the pressure in PWN 
bubbles, the estimated conversion radii fall by a factor of 
30. From Table [TJ we see that the precursor waves would 
in this case be launched in the unstable zone R < 100. 
It seems reasonable to suppose that the shock structure 
will be strongly influenced by such a change. However, 
the possible implications for particle acceleration and the 
resulting non-thermal photon emission remain unknown. 

The inner zone of stability, on the other hand, may 
be of importance for pulsar winds in high pressure envi- 
ronments, such as outflows from a companion star. Sev- 
eral gamma-r ay binaries are likely to fall into this class 
pubusl [20061) but so far only one of them, B1259— 63, 
harbors a pulsar with a measured period. Its spin-down 
power L = 8 x 10 35 crgs _1 , implies ol = 3 x 10 9 , an 
order of magnitude lower than that of the Crab. Assum- 
ing young and middle-aged pulsars hav e a pair multiplic- 
ity k = ol/(4/u) that is at least 10 5 (|Bucciantini et all 
[20TlT) . we find, for B1259-63, /i > 7.6 x 10 3 . The con- 
fining pressure at periastron, and, hence the pressure at 
the conversion point, can be estimated as P conv > 2.5 x 
10~ 4 dyn cm" 2 pubusl[200l . This implies R conv « 1, 
i.e., the conversion point lies roughly at the critical ra- 
dius inside which superluminal waves cannot propagate. 
Since the confining pressure varies over the binary orbit, 
it appears possible that the pulsar wind in this object will 
terminate at r < r c close to periastron, but at r > r c at 
larger binary separations. At the transition points, the 
conversion point passes through the inner stable zone, 
and, according to Fig [3J a stable precursor wave is pos- 
sible for a brief interval of phase. 

5. CONCLUSIONS 

Our main results are the derivation for strong super- 
luminal waves of an analogue of the hydrodynamical 
entropy equation, the demonstration that two kinds of 
wave, confined and freely expanding, propagate in spheri- 
cal geometry, and the identification of two zones of stabil- 



8 



ity of the confined mode in pulsar winds. Strictly speak- 
ing, the latter two results apply only to circularly polar- 
ized modes, but the similarity of the dispersion curves 
for linear and circular modes suggest this is not a serious 
restriction. 

The application of these results to PWN around iso- 
lated pulsars and around pulsars in gamma-ray binaries 
suggests different termination shock structures depend- 



ing on the confining pressure. However, the implications 
for particle acceleration and, hence possible observational 
signatures remain a topic for future work. 

We thank Ioanna Arka, Jerome Petri and Takanobu 
Amano for fruitful discussions. 



APPENDIX 

A. RADIAL PROPAGATION OF LARGE- AMPLITUDE WAVES 

To solve a couple d system of two-flui d an d Maxwell equations in s pherical geometry, we use a standard perturbation 
technique, following I Asseo et al.l (|1984D and iKirk fc Mochol ()2011al lbT). We assume two timescales, on which the wave 
properties change: a fast scale - determined by the pulsar rotation period, and a slow scale, on which the wave 
quantities evolve slowly due to spherical expansion. The "fast" variable is the WKB-like phase of a wave 

*-»('- • <ai) 

where /3 = 1//3 W is the superluminal phase velocity of a wave, /3 W is the subluminal group speed. The "slow" radial 
coordinate is defined by the light cylinder distance tl as 

P = e- , (A2) 

where p ~ 1 and e ~ r-^/r <C 1 is a small parameter. The time and space derivatives, expressed in terms of the new 
variables <f> and p, take the following form: 

d d d uj d lj(3 w d d uj d d 
dt d<j> ' dr c dp c d<j> ' ^ dt € c ^ dp dip ' 

where A = 7— (3 w p\\. Ex pressions (|A3[) are substituted into the equations that govern wave propagation (e.g.. lClemmowl 
Il974bt iKennel fc Pellatlll976| ). and, in addition, all the dependent variables in these equations are expanded in e to 
the first order, i.e., 7 = 7^°- ) + e^ 1 ^ etc. The proper density of each plasma species and the electromagnetic fields 
are expressed in a dimensionless form: n — > nmui 2 /8ire 2 , E — > eE/mcuj, B —} eB/mcui, where E = E y + iE z and 
B = By + iB z are complex quantities (so are the transverse momenta p± = p y + ip z ). 
In the lowest order in e we obtain the equations describing plane waves. Those are the continuity equation: 

A( n (°)A(°))=0 , (A4) 



Faraday's and Ampere's laws: 



&B(°) , k s 

^^ + ^ = °' (A5) 

-/^ + ^W°> P f =0, (A6) 



and momentum/energy conservation multiplied by n: 



d (°) 

n (o) A (e»_*|_ + n (0) Im ( p W S (o)*) = o , (A7) 



o (0) 

n (o) A <o)£P±_ _ n (0) ^(0)^(0) + ip (°) B (o)) = o } (A8) 

n (o) A (o)^ _ n ( 0)Rc ( p (o) E (o)*\ = _ (A9) 
dq> \ J 

Equations (|A5|) and (|A6|) imply B^ = i^ w E^ and n^p^ = —(1 — 0^)dE^ / dcj). Using these relations in equations 
(|A7|) and (|A9[) . we arrive at the result that the quantity 6^ = p^ — /3 W 7 < ' ' ) is phase independent. 
In the first order in e the equation of continuity reads: 

A ( n (0) A (D + „(D A «») + (pV^) = , (AlO) 
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Faraday's and Ampere's laws are: 



-p w — 1- 1—^-, h I 



1 d 



do 



and momentum/energy equations, after multiplying by n, give: 



+ „,, A <o,^L + ln „ )]m + „<.U»>^_ . „<», 

V J o<p do dp 

„ fi 

( n (i) A (o) +n (o) A (i) 



(0) 



(0) 



x ,% (0) 

dA °p±- 



dp 
dp 

(0) 



lpi u 



(0)|2 



n (o) A (o)_^_ _ r n 1 E + i o)i (i) + n ( ) (o) 



D («V°) 



( n (x) A (0) +n (0) A (D) ^ + n^A(°)^ - [^(p^l^+n^^ 



= 



(All) 
(A12) 

(A13) 

(A14) 
(A15) 



To ensure that the first-order quantities have regular (nonsecular) behavior, we impose the condition that they are 
periodic in </>. This suffic es to determine the slow radial dependence of the lowest-order, phase-averaged terms. For 
example, integrating Eq. (|A10[) over </>, we obtain 



1 d 

p 2 dp 







(A16) 



Next, we express the transverse momentum in terms of the fields using Faraday's and Ampere's laws, (|A11[) and (IA12I) . 
which allows us to calculate the phase-averaged forces: 



[nRc(p ± E*)} {1) 
nIm(p ± B*)} {1) 



(A17) 
(A18) 



Substituting those into the equations of motion (|A15[) and (|A13[) , we obtain the conservation of the total energy flux 
and evolution of the radial momentum flux, respectively: 



p 2 dp 



^/ n (oyo) 7 (o )+/3 . 



1 d 
p 2 dp 



P 2 {n^ P f + 1(1 + 01) 



E (0) 




=0 , 


(A19) 


E (0) 


2 >: 




(A20) 



B. CONSERVATION OF (7} 

The analogue of the hydrodynamic entropy equation is obtained by multiplying (|A15[) by /3 W and subtracting it from 
jAH|: 



(>) A (0) +n (0) A (D) ^l +n (o) A (o ) d^ +(3w [nRe{p±E *)f 
V / oq> do 



dd> 



[nlm(p ± B*) 



(i) 



dv {0) 



+ n p w dp Pwn p w dp n p 



(0)|2 

- = . (Bl) 



It is shown in Appendix [A] that the lowest order equations imply that S^ ' is p hase independent, and the first order 
quantities are required to be periodic in 4>. Therefore, phase averaging of (|B1[) immediately cancels first two terms, 
and the remaining part takes the form 



+ (»' M,° ) ^)-^H," , ^)-;<" l ° > ^ , l 2 



, (B2) 
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where we have used the expressions for phase-averaged forces (|A17[) and (|A18[) . In the following we omit superscript 
"0" and further simplify the equation (|B2[) by the momentum relation p 2 = 7 2 — p\ — 1 , so that it reads 

m + _l s Jfi +m ^ M>= . (B3) 

llP 2 7 2 dp dp 2\ dp / p 

The first term in (|B3[) cancels the last term, since: 

< e2) - r ^ = c r ^% = r - * ■ < b4 > 

(0), 



Here in the second step we use the equation (|A8[) and (|A5|) . which together imply i?( ) = dp\_ /d<p, and we integrate 
once by parts. Further simplification follows from Eq 
prove that the second term in the entropy equation (|B3 



6|): d 2 p±/d(j) 2 = dE^°'/dc/) = —j 2 / np±. The same relations 
cancels the fourth term: 



d P 2 ±\ f^^f.^ dP*± *dP±\ 1 [*",,( BE dp* x dE* dp A 



n > = / d(j)\ npi_——^- + np*L 



dp I Jo \ dp dp J J \ d<p dp dcj) dp 



Thus, Eq. (|B3[) implies that independently of the wave polarization the phase-averaged Lorentz factor of the particles 
stays constant during the radial expansion: 

^=0. (B6) 
op 

C. LINEAR. POLARIZATION 

For linearly polarized waves all the quantities are phase dependent. It is convenient to use the phase variable 
y = E/Eq, chosen to be y = 1 for the phase <p, for which the electric field takes the maximum value Eq. The nonlinear 
dispersion relation follows from the periodicity requirement 

(CI) 



(C2) 



t Jo \dy/d<f>\ 

and Ampere's law takes the form 

dy = ^p7wP-L A o^po7w P± 
d(j> ~ E oj 2 _ E uj 2 A ' 

where cu 2 = 87re 2 no/m, and the continuity equation (jA4[) implies that nA = noAo is phase-independent. In radial 
evolution equations all the phase averages have to be calculated explicitly, i.e., 

n-tt^mk ' (C3) 

Dependence of plane-wave particle momenta 7, p±, p» on y one obtains by integrating equations of motion, e.g., 
IKennel k PeTIatl ((1971 : 

7 = 7o+ 9 ff° 2 2 (l~V 2 ) , (04) 

Ul 2 E 2 j3 w . 2 \ /V^cl 

p\\ = p\\o + 0l ,2 a 72 ( x - y ) . ( C5 ) 

P± = (7 2 -Pf-I) 1/2 , (C6) 

where 7q = p| + 1 . The set of equations which describes radial propagation of a linearly polarized wave with variables 
no, P\\o, Eq, p±, 70, f3 w is given by equations ((SJ), ©, and (JUJ), in which the phase averages have to be calculated 
explicitly according to Eq. (|C3|) . together with (|C1[) and (|C6|) . 



D. THE OUTER STABILITY ZONE 

The Hugoniot curves are determined by Eqs. (|17[) and ((2"2"j) , where v w {R) — v is a constant, related to /i and tr by 
fU). These equations are complemented by (fl9|) and momentum relation 7 = (1 + p\ + p 2 ) 1 ^ 2 - Finding p± from ([22]) 
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and substituting to (|17|) . we obtain a biquadratic equation for radius R as a function of the wave Lorentz factor y w . 
For the confined mode branch the larger root is relevant. The largest term in the expansion in e ~ /1 -1 / 3 ~ er -1 / 2 
reads: 



I?" 



7 2 -2 7 2 (2 7w -l) ( 7 l-l) 1/2 




h- 1 


-2 7 3 (27 2 -l) ( 7 w" V / 7 2 -l) 


4 7 w (2 7 w - ^Jll - l) 


-87i( 


7w - \/7w ~ ! 


r— 1 



(Dl) 



The lower bound on 7 W , and thus R, for which the Hugoniot curve implies stable initial conditions for a wave at 
launch, is obtained from the stability condition ([T]) , in which the particle momenta are transformed from the H-frame 



to the lab frame: 7' = 7 W (7 — /?wP||), pj| = Jw(p\\ — /?w7)- Expanding in e ~ fi 
radii R ~ e _1 , in the lowest order Eq. ^ takes the form 



-1/3 



2(-2 7 3 + 7 4) 



> 



This implies the lower bound 7 W ss 2, and thus Eq 
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